Statistical Approach to the Cosmological-Constant Problem on Brane Worlds 
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We investigate the physically accepted solutions of general Braneworld scenarios, scanning uni- 
formly the associated parameter space. Without making any further assumptions we find that 
solutions which give "small" Hubble parameters on the physical brane, and therefore "small" ef- 
fective cosmological constants on the 4D Universe, are far more probable than those with "large" 
ones. Eventually, their distribution tends to the 5-function in the limit of continuous covering of 
the parameter space. 
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Why is the cosmological constant of the observable 
Universe many orders of magnitude smaller than ex- 
pected within the context of ordinary gravity and quan- 
tum field theory? Since Einstein (1917) many authors 
have tried to give a plausible and coherent answer 
and references therein) . Most of the earlier attempts were 
within the general relativity context, while other authors 
were led to the obscure anthropic principle 2] . An outlet 
of different origin arose since Zumino showed that unbro- 
ken supersymmetry leads to vanishing vacuum energy Q . 
Unfortunately, and despite many attempts in the context 
of supergravity (or non-conventional gravity) 0, 0] and 
superstrings [61 , in the real world we know of no mecha- 
nism that can produce an effective cosmological constant 
within the remarkable small experimental limits, and es- 
pecially to keep it unaffected by quantum fluctuations. 
Another interesting approach was followed by Hawking 
and later by Coleman Q, showing that a distribution 
of values for the effective cosmological constant with an 
enormous peak at zero could arise. Although the whole 
approach raises many doubts [H, the idea of restricting 
the cosmological constant to small values through statis- 
tical considerations is interesting and consists the basis 
of the present work. Finally, the approaches to the prob- 
lem which remain in the conventional 4D framework suf- 
fer from loss of generality [8( , and the use of cosmological 
observations on behalf of small values pays the price of 
adding various non-general assumptions [9j. 

The discussion was reheated after the arrival of obser- 
vational data indicating that there is a positive cosmo- 
logical constant in the universe [loj ]. At the same time 
the Braneworld scenario appeared in the literature 11] . 
Amongst other consequences, it offered new paths of pro- 
ducing a zero or small (effective) cosmological constant 
on the physical brane consisting our Universe (13 , Il3j |. 
but none could offer a robust and undoubted justification. 
Most of these works re- introduced in the new Braneworld 
context (of various versions and dimensions) the old fine- 
tuning procedure of conventional cosmology. 



In conclusion, there is not a convincing and general so- 
lution to the cosmological constant problem. Definitely, a 
complete apprehension will arise from a fundamental the- 
ory of nature and the breaking of its symmetries. Since 
this theory is unknown at present, effective or even phe- 
nomenological descriptions are probably preferable, al- 
though incomplete. The minimum requirement is that 
of remaining general and fine-tuning free. In this Letter 
we follow this direction in the context of Braneworld sce- 
nario, and we explore the physically accepted solutions 
in the parameter space, without making any additional 
assumptions. Obviously this approach cannot act as a 
solution to the cosmological constant problem. However, 
it can be enlightening and provide some arguments for 
its possible elaboration in statistical terms. 

We consider a general class of Braneworld models, 
characterized by the action: 
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-773- is a 5D gravitational constant, and all 
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quantities are measured in units of A/5. The first term 
describes gravity in the five dimensional bulk space, the 
second term corresponds to a minimally coupled bulk 
scalar field with the potential V{<j>) and the last term 
corresponds to two 3 + 1 dimensional branes, which con- 
stitutes the boundary of the 5D space. We include a po- 
tential term U (0) for the scalar field at each of the two 
branes, denoting with 7 the associated induced metric 
and with K its extrinsic curvature. The square brackets 
denote the jump of any quantity across a brane ([Q] = 
Q(y+) — Q(y-)). Assuming S 1 /^ symmetry across each 
brane we restrict our interest only in the interbrane space 
(HQlll, (we use [Q'] = 2Q'(0+), [Q']t = -2Q'(1")). 
As usual the two branes are taken parallel, y denotes the 
coordinate transverse to them and isometry along three 
dimensional x slices including the branes, is assumed. 
For the metric we choose the conformal gauge (lit ]: 

dg 2 = e 2B(t,y) f_ dt 2 + dy 2^ + e 2A(t,„) ^2 _ (2) 

This metric choice, along with the residual gauge free- 
dom (t, y) — > (t',y') which preserves its form, allows 
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us to "fix" the positions of the branes. Without loss 
of generality we can locate them at y = 0,1, having in 
mind that their physical distance at a specific time is 
given by D(t) = dye B ^ 0. The non-trivial five- 
dimensional Einstein equations consist of three dynami- 
cal: 



A-A" + 3i 2 - 3A' 2 = -e 2B V 



B - B 



3 A 2 + 3A' 2 
I -<t>" + SA<}> - 3A'4>' 



2L + - - e 2B y 
2 2 3 

e 2B V^ = 0, 



(3) 



and two constraint equations: 



A' A + B'A + A'B - A' 



1 ■ 



2A' 2 - A'B' + A" - A 2 -AB 



6 6 



e 2B V\A) 



where primes and dots denote derivatives with respect 
to y and t respectively. Additionally, from the boundary 
terms in the action for the branes we obtain the following 
junction (Israel) conditions: 

[A']=t\uc b : [B']= T \ue B 1 m = ±e B U#, (5) 

where the upper and lower signs refer to the branes at 
y = 0, 1, respectively. For the bulk and brane potentials 
we assume: 

V(<f>) = im^+A, U l {cj>) = \M*{cj>i-cj l ) 2 + A,. (6) 

This is a general form for VUf>), consistent with the brane 
stabilization mechanism [16j |. A being the 5D bulk cos- 
mological constant. For Ui{<f) the quadratic form is also 
general 15|, incorporating the brane tensions A^ as well 
as the value 4>i of <f> on the i-th brane. Finally, the in- 
duced AD metrics of the two ("fixed" -position) branes 
in the conformal gauge are simply given by: ds 2 — 
— dr 2 + a 2 (r) dx 2 , with Tj = e Bi t and at = e Ai the proper 
times and scale factors of the two branes [17[ ■ Thus, for 
the Hubble parameter on the branes we acquire 



H, 



1 da 
a dr 



~ B 'A„ 



(7) 



and we identify the physical brane with the one at y = 0. 

The described model is quite general and includes the 
full spacetime evolution of the 5D Braneworld. In order 
to investigate some of its consequences we first examine 
the stationary sub-class of solutions, assuming that: 



B{t,y) -> B(y), A(t, y) — > B(y) + Ht. 



(8) 



This choice leads to a fixed bulk geometry (<j> is time 
independent), and maximally symmetric (de Sitter or 



Minkowski) branes: ds 



(dy 2 -dt 2 +e 2Ht d X 2 ), 



where the Hubble parameter on the physical brane is sim- 
ply Hq = e~ B (°)H. If we set V{4>) = we can acquire 
analytical solutions for the system of equations © and 
@, depending on the sign of H 2 . For H 2 > and setting 
H = we get: 



B(y)=B(0) + ±log 



B'(Q) 
H 



sinh 3Hy + cosh 3Hy 



0(l/) = 0(O)--^=loE 



,3Hy 



u) (1 + u) 



( e 3Hy +U }(1_ u ) 



(9) 



where 



§r^§, and always \B'(0)\ > H for the 

existence of a solution. When H 2 < and setting 
the solutions can be obtained from (J9j) 
I18ll. We mention that 



9 = 

under the replacement H —> it 
we do not make additional assumptions, excluding this 
case (H 2 < 0) by hand, since we investigate any possi- 
ble solution. Finally, expressions ([9]) have to satisfy the 
boundary conditions ([5]), and eventually .0(0), -0'(O) and 
4>(0) and H are given in terms of the six parameters of 
the model (Mo, Ao, <jq, Mi, Ai, ui). 

In order for the solutions to be physically accepted, we 
have to assure that no "naked" singularities are present 
in the physically meaningful interbrane space, i.e between 
y = and y = 1. As it is shown in [l9| . the singulari- 
ties of the aforementioned solutions are indeed true (and 
not just coordinate ones) since they possess incomplete 
geodesies reaching a point of divergent curvature. An 
equivalent way to reveal the true nature of these singular- 
ities, is by calculating explicitly any of the 5D curvature 
invariants such are the Ricci scalar, the Kretschmann 
scalar and the Weyl tensor, and showing that they di- 
verge. This has been performed in 0, 0, [2(|, and has 
been straightforwardly verified by us that these curvature 
invariants do diverge in the singularities of the solutions 
(J9j) , offering a robust proof of their true (naked) nature. 

For H 2 < case, the requirement of singularity ab- 
sence forces the physically accepted solutions to satisfy: 



39 — 7r < arctan 



B'(0) 



< 0. 



(10) 



Thus, this relation provides a narrow and absolute win- 
dow for 9, excluding areas of the 6D parameter space 
leading to 9- values larger than 7r/3. For the H 2 > case 
the avoidance of singularities leads to: 



H < \B'{0)\ < Hcoth(3H), 



(11) 



and since B'(0) and H are functions of the six parameters 
of the model, the window (|11[) gives the limits of the 
sub-space of the 6D parameter space that corresponds 
to accepted solutions. Inequality (fTTj) becomes actually 
an equality for large H values (H > 2 in the units we 
use) since both edges practically coincide. The resulting 
relation acts as a constraint in the 61? parameter space, 
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and therefore the parameter sub-space which corresponds 
to physically accepted solutions with large H is of lower 
dimensionality, i.e it has been reduced to a sub-surface: 



3x10 5 n 



H = \B'(0)\ = # coth(3tf). 



(12) 



In conclusion, from the above analytical calculation it is 
implied that the demand for absence of divergences in the 
solutions restricts the probable H 2 values (and therefore 
Hq ones) to a narrow window around zero (solutions with 
large H have zero measure). The additional requirement 
of fulfilling the junction conditions on the second brane 
makes the parameter sub-space, leading to a solution, 
in-homogenous and non-compact 
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Let us verify numerically the acquired results of the 
stationary case with zero bulk potential. The main dif- 
ficulty of solving the equation system (J3|)-([6|) is that one 
has to satisfy the constraints (J4j) and boundary conditions 
([5]) at t = 0, ensuring that no divergencies are present in 
the interbrane space, task which proves to be extremely 
hard in general. The restriction to t — is easily justi- 
fied since the dynamical equations ([3j> preserve the con- 
straints f4|, and furthermore for the stationary case at 
hand the boundary conditions J5]| are time-independent. 
The method we use is the following: We first choose ran- 
domly the values of the 6 model parameters, uniformly 
distributed in a 6D hyper-cube. The obtained results 
do not depend on the hyper-cube's size, but on its effec- 
tual covering (number of parameter multiplets used in the 
calculation). We use the globally convergent Schmelcher- 
Diakonos algorithm [Hj to solve the transcendental equa- 
tion system (J4j) — (JSJ) , with accuracy 10~ 13 . Thus, all ex- 
isting solutions are easily obtained, which is not the case 
if one uses Newton's method [Hj], where a good initial 
guess, of a (relatively small) interval containing the root 
to be found, is needed as an input. We find hexads of pa- 
rameters corresponding to acceptable solutions, and cal- 
culate Hq through Hq = e~ B ^H. In fig. QJ, we present 
the Hq histogram for N = 10 6 accepted solutions, and 
we observe that p(Hq) is strongly restricted around zero 
according to windows (p~0|) and (fl~Tj) . A numerical fit spec- 
ifies that for Hq > the aforementioned histogram is an 
exponential (coefficient of determination R 2 « 0.999): 

p(H 2 )~e- aH K (13) 

with exponent a ~ 3. In order to provide an additional 
verification, in fig. [TJd we depict H versus -B'(O) for so- 
lutions that have H 2 > 0. We observe that for small H 
(H < 0.5) inequality (fTT|) holds, and it gradually trans- 
forms to equality (|12p as H grows, i.e the parameter sub- 
space that corresponds to accepted solutions becomes 
tinner and thiner and finally results in a sub-surface. 

We mention here that in our analysis we follow the 
"anti-fine-tuning" approach to the problem, determining 
first, unbiased (i.e randomly), the six model parameters 
and then investigating the solution features. This is op- 
posite to the usual method in the literature, which is to 
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FIG. 1: (Color online) In graph a) we present the Hq his- 
togram for the accepted solutions. In graph b) we depict H vs 
-B'(O) for solutions with H 2 > 0. Both graphs correspond to 
the stationary case with zero bulk potential, and each edge of 
the parameter space hyper-cube extends from -10 3 to 10 3 . 



fulfil the boundary conditions in the first brane at y = 0, 
choose an arbitrary H, solve the constraints ([4J towards 
y = 1 and then fit the brane parameters of the second 
brane so that the corresponding junction conditions are 
fulfilled. Definitely this fine-tuning procedure has the 
disadvantage of neglecting several solutions, allowing for 
an apparent uniform distribution of the H value, con- 
trary to the fact that the considered solution space is a 
sub-surface of the full 6D parameter space. The source 
of the misleading behavior is the over-determined nature 
of the problem, i.e the existence of more relations and re- 
quirements to be fulfilled (including conditions that are 
not expressible through a strict equation, such is the ab- 
sence of naked singularities) than the number of model 
parameters. In such cases our method is better in order 
to reveal the full characteristics of the parameter space 
and its sub-space that corresponds to solutions. 

After this verification of the analytical calculations we 
can proceed to numerical solutions for the stationary case 
with full bulk potential. Again, we randomly choose the 
values of the 8 model parameters (Mo, Ao, <Jq, M\, Ai, g\, 
m and A), from a uniform distribution, and we solve the 
transcendental equation system dl}-©. We mention here 
that we do not use any a priori restrictions on the param- 
eter values arising from physical arguments (for example 
A < in order to acquire a 5D AdS bulk as required by 
many authors), since we want to be as general as possible. 
In fig. [5h. we present p{H^) for N = 10 6 such solutions, 
and a numerical fit for Hq > gives again an exponen- 
tial {R 2 « 0.998) of the form $T5§ with exponent a » 2.9. 
Although an analytical treatment is not possible in this 
case, the qualitative results are similar to those obtained 
for zero bulk potential. Indeed, it seems that "small" Hq 
arise from a (non-homogenous and non-compact) sub- 
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an explicit relation for A e fi in terms of the model param- 
eters and then force this expression to be small. On the 
contrary we have just searched for general parameter oc- 
tads that correspond to solutions and for these octads we 
have calculated Hq. In fig. [3] we present the histograms 
of the model parameters that correspond to N = 10 6 ac- 
cepted solutions, for the stationary case with full bulk 
potential. As we can clearly see, apart from a small ten- 
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FIG. 2: (Color online) In graph a) we present the Hq his- 
togram for the accepted solutions. In graph b) we show the 
histogram p(Hq) of Hq, for four different numbers of physi- 
cally accepted solutions, where the corresponding curves have 
been normalized to area 1. The dotted lines correspond to lin- 
ear fits. In graph c) we depict the dependence of exponent a, 
of the exponential form H13V . on the physically accepted so- 
lution number N. Numerical extrapolation gives a — > oo for 
N — > oo. All graphs refer to the stationary case with full bulk 
potential. 



space of the parameter space, which becomes thinner for 
"larger" Hq, resulting eventually to a sub-surface. 

The natural arising question is what factors affect the 
exponent a of relation (|13p . We can show numerically 
that it is determined by the statistics and increases with 
the number of solutions N. In fig. [2]d we depict p(Hq) vs 
Hq, for four different-TV cases. We clearly see that the 
exponential form decays faster for larger N. In fig. [2J; 
we present the values of the exponent a, arising from 
numerical fits of the form of (Ti~3")) , for various N values. 
Indeed, a increases with N, and the increasing rate is ro- 
bust although quite slow. Numerical extrapolation gives 
a — > oo for N — > oo. This is a favoring result since it 
implies that covering continuously and uniformly the pa- 
rameter space, the physically accepted solutions compose 
for p(Hq) the 5- function 5(Hq). Solutions with large Hq 
are not impossible, bur arise from a sub-surface of smaller 
dimensionality. In the continuum limit, large Hq solu- 
tions are of zero measure while the ones with Hq 
are by far the most probable (note that although we are 
using Ms units, the presence of 5(Hq) leads to Hq ps in 
every units). As it is well known, Hq in our Universe is 
related to the (effective in a 5D model) 4D cosmological 
constant of our world. Since |A e ff| < \Hq\ we finally ac- 
quire p(|A e ff|) ~ <5(|A c ff|), compatible with observations. 

The distribution for A ff obtained above is free of both 
explicit and implicit fine-tuning, since we have not used 
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FIG. 3: (Color online) The histograms of the model param- 
eters m, A, Ao, Ai, Mo, Mi, cro, o\, that correspond to 
N — 10 6 accepted solutions, for the stationary case with full 
bulk potential. 



dency of m to small values, there is no specific preference 
to particular values. To quantify this behavior we calcu- 
late the correlation dimension of the sub-space of the 
parameter space that corresponds to accepted solutions. 
Denoting by Fy the 8D Euclidean distances between two 
solution-octads i and j: 

\nj\ 2 = (mi - mjf + (A, - Kjf + (A 0l - Xqj) 2 + 
+(A M - Xijf + (M w - Mq 3 ) 2 + (M Xi - M Xi ) 2 + 

+(<j i - vo 3 ) 2 + (<rii - &i 3 ) 2 , (14) 

we use the standard approach introduced by Grassberger 
and Procaccia [22| to calculate the correlation integral 



C(r) 



N(N - 1) 



$>(r 



(15) 



in order to obtain the correlation dimension (j. In fig. 0] 
we depict C(r) versus r, in the region of intermediate r 
where the statistics allows for the power-law form of re- 
lation (p~5|) . We observe that it clearly follows (|15p with 
exponent (3 = 7.2 ± 0.3. Such a value (within statistical 
errors) is very close to the value (3q — 8 expected for a 
complete covering of the 8D space. To estimate the effect 
of the finite statistics we have repeated the calculations 
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FIG. 4: (Color online) C(r) versus r, in the region of interme- 
diate r 's. It follows the power-law form of \15\) . with exponent 
(3 = 7.2 ±0.3. 



using the same number of random octads within the cor- 
responding hypercube (i.e without taking into account 
any constraints induced by the requirement of physically 
accepted solutions), and we have found a typical value 
for the correlation dimension (3 sa 7.4 ± 0.2. Therefore, 
our results should be compared with this value and not 
with the value (3q = 8. Thus, based on this fact, we 
safely conclude that the obtained iJQ-histogram and A e ff- 
distribution, are fine-tuning-free. 

In the above analysis we have been restricted to sta- 
tionary solutions of the form ([5]) without bothering about 
time evolution, since in this case Hq — Hq(0) — const. 
by construction. The question is what can be said about 
the full dynamics of equations (O and ((3]), where Hq can 
be varying. Fortunately, the solutions of the full dynam- 
ics seem to consist of stationary or slowly varying ones 
14, [l5j]. In such an evolution H (t) can change by at 
most one order of magnitude. This behavior justifies the 
above stationary treatment, since the obtained results for 
t = will not change qualitatively for t > 0. 

In order to present this behavior in a more transparent 
way we investigate numerically the full dynamical system 
(|3|), ([4} and ((U, that is without imposing the ansatz ((8]). 
In particular, we first perform the procedure described 
above for t = (or equivalently for initial brane proper- 
time r = 0), in order to extract N = 5 x 10 4 parameter 
octads that correspond to solutions, and for these octads 
we calculate H o (0) and p(i?o(0)). For these parameter 
octads we evolve the system according to <j3j) and ([5]), 
and thus for each brane proper-time we obtain the cor- 
responding histogram p(_ffg(r)). We mention that for 
a significant portion of the initially acquired solutions 
(~ 20%), the subsequent time evolution leads to diver- 
gences, since the singularities transit from outside to in- 
side the inter brane space. This behavior was also found 
in 



15j . Thus these parameter octads are neglected. 



In fig. Owe present the initial histogram p (iJp (0)) (up- 
per graph) and the histogram p(^H^(r)^ for r = 10 3 in 
M5 units (lower graph). As we observe, indeed the vari- 
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FIG. 5: (Color online) The initial histogram p (Hq(0)) (upper 
graph) and the histogram p (Hq(t)) for r — 10 3 in M5 units 
(lower graph), for N = 5 x 10 4 physically accepted solutions 
of the full dynamical system O) and with full brane 
and bulk potentials. 



ation of Hq(t) from its initial value Hq(0) is less than an 
order of magnitude for the significant portion of cases. 
Note also that this variation can be either an increase 
or a decrease. Thus, in summary, the final histogram is 
only slightly wider. To quantify this deformation, in fig. [6] 
we present the aforementioned histograms in logarithmic 
scale, normalized to area 1. As we can see, time evolu- 
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FIG. 6: (Color online) The initial histogram p (Hq (0)) and the 
histogram p (Hq(t)) for r — 10 3 in Ms units, for N = 5 x 
10 4 physically accepted solutions of the full dynamical system 
(00), {^P an d EP; with full brane and bulk potentials. The 
histograms have been normalized to area 1 and the dotted lines 
correspond to linear fits. 



tion does not change our results significantly. Thus, in 
the limit N — > 00 the <5-function-form p(|A ff|) ~ 5(|A e ff|) 
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will be preserved. 

One enlightening interpretation of the aforementioned 
result is obtained by a parallelism with the equilibration 
mechanism of an isolated macroscopic system through 
entropy maximization. "Peculiar" states are not impos- 
sible, but there is extremely small probability for their 
realization. In a similar way, in braneworld models, the 
large majority (in the continuum limit "almost" all) of 
the parametric combinations that correspond to physi- 
cally accepted solutions, lead to an extremely small cos- 
mological constant. 

In conclusion, using a statistical approach, we have 
shown that in the context of a general Braneworld sce- 
nario the acceptable solutions of the Einstein equations 
are compatible with a (^-function distribution for the AD 
effective cosmological constant. Additionally, within this 
treatment, A e s is small and remains so even if the model 
parameters vary by many orders of magnitudes due to 
successive cosmological phase transitions. Of course the 
above analysis cannot at all act as a solution to the cos- 
mological constant problem. It only provides some argu- 
ment for its possible elaboration in statistical terms. 
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